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Abstract 

We present a simple technique for avoiding physically spurious eigenmodes that 
often occur in the solution of hydrodynamic stability problems by the Chebyshev 
1 collocation method. The method is demonstrated on the solution of the Orr- 

O-i Sommerfeld equation for plane Poiseuille flow. Following the standard approach, 

O • the original fourth order differential equation is factorised into two second-order 

equations using a vorticity-type auxiliary variable with unknown boundary val- 
ues which are then eliminated by a capacitance matrix approach. However the 
elimination is constrained by the conservation of the structure of matrix eigen- 
value problem, it can be done in two basically different ways. A straightforward 
application of the method results in a couple of physically spurious eigenval- 
c/2 ■ ues which are either huge or close to zero depending on the way the vorticity 

boundary conditions are eliminated. The zero eigenvalues can be shifted to 
any prescribed value and thus removed by a slight modification of the second 
approach. 



1. Introduction 

Spectral methods are known to achieve exponential convergence rate 
which makes them particularly useful for solving numerically demanding dif- 
ferential eigenvalue problems which arise in hydrodynamic stability analysis 
[l^ . Unfortunately, besides providing accurate and efficient solutions for a cer- 
tain number of leading eigenvalues, spectral methods often produce physically 
spurious unstable modes, which cannot be removed by increasing the numerical 
resolution [f| . For detailed discussion of these modes we refer to Boyd [2] . Such 
physically spurious eigenvalues can appear in all types of spectral methods in- 
cluding Galerkin [16j, tau [J] and collocation approximations [3|, unless some 
kind of ad hoc approach is applied to avoid them. In the Galerkin method, 
spurious eigenvalues can be removed by using the basis functions also as the 
test functions instead of separate Chebyshev polynomials [r7| . A number of ap- 
proaches avoiding spurious eigenvalues have also been found for the tau method 
0, [Hi [HJ • The same can be achieved also for the collocation (or pseudospec- 
tral) method by using two distinct interpolating polynomials [9]. Following the 
approach of McFadden et al. [nj for the tau method, Huang and Sloan [9] 
use a Lagrange interpolating polynomial for second-order terms which is by two 
orders lower than the Hermite interpolant used for other terms. The choice of 
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the latter polynomial depends on the particular combination of the boundary 
conditions for the problem to be solved [H p. 493]. 

The objective of this paper is to present a simple method avoiding spurious 
eigenmodes in the Chebyshev collocations method which uses only the Lag- 
range interpolating polynomial applicable to general boundary conditions. Our 
approach is based on the capacitance matrix technique which is used to elim- 
inate fictitious boundary conditions for a vorticity-type auxiliary variable. The 
elimination can be performed in two basically different ways which respectively 
produce a pair of infinite and zero spurious eigenvalues. The latter can be shifted 
to any prescribed value by a simple modification of the second approach. The 
main advantage of our method is not only its simplicity but also applicability 
to more general problems with complicated boundary conditions. 

The paper is organised as follows. In the next section we introduce the Orr- 
Sommerfeld problem for plane Poiseuille flow, which is a standard test case for 
this type of method. Section [3] presents the basics of the Chebyshev collocation 
method that we use. The elimination of the vorticity boundary conditions, 
which constitutes the basis of our method, is performed in Sec. 0] Section [3] 
contains numerical results for the Orr-Sommerfeld problem of plane Poiseuille 
flow. The paper is concluded by a summary of results in Sec. [5J 

2. Hydrodynamic stability problem 

The method will be developed by considering the standard hydrodynamic 
stability problem of plane Poiseuille flow of an incompressible liquid with density 
p and kinematic viscosity v driven by a constant pressure gradient Vpo = —&xPo 
in the gap between two parallel walls located z — ±h in the Cartesian system 
of coordinates with the x and z axes directed streamwise and transverse to the 
walls, respectively. The velocity distribution v(r,t) is governed by the Navier- 
Stokes equation 

d t v + (v ■ V)v = -p~ 1 Vp + vV 2 v (1) 

and subject to the incompressiblity constraint V • v = 0. Subsequently, all 
variables are non-dimensionalised by using h and h 2 jv as the length and time 
scales, respectively. Note that instead of the commonly used maximum flow 
velocity, we employ the viscous diffusion speed v jh as the characteristic velocity. 
This non-standard choice will allow us to test our numerical method against the 
analytical eigenvalue solution for a quiescent liquid. 

The problem above admits a rectilinear base flow Vq(z) = Reu(z)e x , where 
u(z) — I — z 2 is the parabolic velocity profile and Re = U$h/v is the Reynolds 
number defined in terms of the maximum flow velocity Uq = 2Poh 2 / pv. Stabil- 
ity of this base flow is analysed with respect to small-amplitude perturbations 
Vi(r,t) by searching the velocity as v = Vq + V\. Since the base flow is invariant 
in both t and x — (x, y), perturbation can be sought as a Fourier mode 

vi(r,t) = v(z)c xt+ik - x +c.c, (2) 

defined by a complex amplitude distribution v(z), temporal growth rate A and 
the wave vector k — (a,/3). The incompressiblity constraint, which takes the 
form D ■ v =0, where D = e z j^ + ik is a spectral counterpart of the nabla 
operator, is satisfied by expressing the component of the velocity perturbation 
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in the direction of the wave vector as u u = e M • v = ik" 1 ^' , where e M = k/k and 
k = \k\. Taking the curl of the linearised counterpart of Eq. (fjQ) to eliminate the 
pressure gradient and then projecting it onto e z x e M , after some transformations 
we obtain the Orr-Sommerfeld equation 

XD 2 w = D 4 w + iaRe(u" - uD 2 )w, (3) 

which is written in a non-standard form corresponding to our choice of the 
characteristic velocity. Note that the Reynolds number appears in this equation 
as a factor at the convective term rather than its reciprocal at the viscous term 
as in the standard form. As a result, the growth rate A differs by a factor Re 
from its standard definition. The same difference, in principle, applies also to 
the velocity perturbation amplitude which, however, is not important as long as 
only the linear stability is concerned. The no-slip and impermeability boundary 
conditions require 

w = w' = at z = ±1. (4) 

Because three control parameters Re and (a, j3) appear in Eq. © as only 
two combinations aRe and a 2 + (3 2 , solutions for oblique modes with (3 ^ 
are equivalent to the transverse ones with /3 = and a larger a and, thus, a 
smaller Re which keep both parameter combinations constant @. Therefore, it 
is sufficient to consider only the transverse perturbations [k = a). 

The first step in avoiding spurious eigenvalues in the discretizied version of 
Eq. ([3]) to be derived in the following section is to represent Eq. ([3]) as a system 
of two second-order equations [8] 

AC = D 2 ( + iaRe(u"w-u(), (5) 
C = D 2 w, (6) 

where C is a vorticity-type auxiliary variable which has no explicit boundary 
conditions. 



3. Chebyshev collocation method 

The problem is solved numerically using a collocation method with iV + 1 
Chebyshev-Gauss-Lobatto nodes 

%i = cos (in/N) , i = 0, ■■■ ,N. (7) 

at which the discretizied solution (w,Q(zi) = (wi,Ci) — ( W 5 C) an d its deriv- 
atives are sought. The latter are expressed in terms of the former by using 
the so-called differentiation matrices, which for the first and second derivatives 
are denoted by D$ and with explicit expressions given in the Appendix. 
Requiring Eqs. (I5I6P to be satisfied at the internal collocation points < i < N 
and the boundary conditions Q at the boundary points i = 0, N, the following 
system of 2N algebraic equations is obtained for the same number of unknowns 

ACo = ACo + Bd+go, (8) 
Co = Aw , (9) 
Oi - Cw , (10) 
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where the subscripts and 1 are used to denote the parts of the solution at the 
inner and boundary collocation points, respectively; wi = Oi due to the first 
boundary condition Q and 

gi = \aRe(u"wi - Ui&). (11) 

The matrices 

Ai d = (D 2 ) i;j , 0<(i,j)<N, (12) 
Bij = (£>%•, 0<i<N,j = 0,N, (13) 

represent the parts of the collocation approximation of the operator 

(D 2 ) id = D<$-c( i I iti (14) 

using the inner and boundary points, respectively; I^j is the unity matrix. 
Equation (jTUJ) is a discretizied version of the second boundary condition Q 
imposed on w' which is defined by the matrix 

C.., />;';• i = 0,iV;0<i<iV. (15) 

Our goal is to reduce Eqs. (fSBT0|) to the standard matrix eigenvalue problem 
for wq. First, £ is eliminated from Eq. (jS]) by using Eq. ([9]). which results in 

AAw = A 2 w + Bd + g . (16) 

Next, we can use Eq. (fTU]) to eliminate d from the equation above. This, as 
shown in the next section, can be done in two basically different ways. 



4. Elimination of the vorticity boundary values 

In order to eliminate d from Eq. (JTHJ) using Eq. (fTU)) we employ a modified 
capacitance (or influence) matrix method. For the basics of this method, see 
[131 p. 178] and references therein. Modifications to the method are due to 
the structure of the matrix eigenvalue problem which needs to be conserved in 
the elimination process. The general capacitance matrix approach suggests to 
express Wo from Eq. P^|) and then to substitute it into Eq. (fTUl) . which then 
would result in a system of linear equations for £ 1 . However, as noted above, the 
elimination procedure must be linear in A for the eigenvalue problem structure 
to be conserved. It means that Wo can be expressed either from the right or 
left hand side of Eq. (|T6|) but not from the combination of both sides as in the 
standard capacitance matrix approach for the time stepping schemes. 

Our first approach is to express Wo from the r.h.s. of Eq. (UHJ) by inverting 
A 2 and then substituting it into the boundary condition (|10j) . which results in 

Cw -CA- 2 (AAw -BC 1 -g )-0 1 . (17) 

Next, solving the equation above for 

C 1 = (CA- 2 B)- 1 CA- 2 (AAw -g ) (18) 

and substituting it into Eq. (|16[) . we obtain 

AEAw = (A 2 + EG)w , (19) 
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where Gwq = g and 

E = I B(CA 2 B) X CA 2 . (20) 

It is important to notice that EB = OB, which means that E is singular. Namely, 
it has a zero eigenvalue of multiplicity two corresponding to two eigenvectors 
represented by the columns of B. Representing Eq. (|19[) as 

(A 2 + EG)- 1 EAw = A -1 Wo, (21) 

which is a standard eigenvalule problem for A -1 , it is obvious that zero eigenval- 
ues of E result in two zero eigenvalues A , which in turn correspond to infinite 
eigenvalues A of the original Eq. (|19p . A way to avoid these spurious eigenvalues 
is described below. 

Alternative approach how to eliminate £ : is to express Awo from the l.h.s. 
of Eq. (|16p by inverting A and then substituting it into the boundary condition 
(|TDj) . which results in 

ACw = CA- 1 (A 2 w + BC 1 +g ) = 1 . (22) 

This equation can be solved for £j similarly to Eq. (fT7| as 

Ci = (CA- 1 B)- 1 CA- 1 (A 2 + G)w , (23) 

which substituted in Eq. (fl"6|) leads to 

AAw = F(A 2 + G)w , (24) 

where the transformation matrix 

F = I B(CA _1 B) *CA 1 (25) 

is singular with two zero eigenvalues because it satisfies FB = OB similarly to 
the E considered above. In contrast to the previous eigenvalue problem defined 
by Eq. (|19p . now the singular transformation matrix appears on the r.h.s. of 
Eq. and thus it produces two zero rather than infinite eigenvalues A. 

It is important to notice that zero eigenvalues represent an alternative solu- 
tion to Eq. (|22[) . which can be satisfied not only by the boundary condition <jTUJ) 
but also by A = 0. Consequently, these spurious eigenvalues can be shifted from 
zero to any value Ao by subtracting AoCwo from both sides of Eq. ([22]) . which 
obviously does not affect the true eigenmodes satisfying Eq. (jTUJ) . As a result 
we obtain 

d = (CA- 1 B)- 1 CA" 1 (A(A - A I) + G)w , (26) 
which substituted in Eq. (fT()]) leads to the following standard eigenvalue problem 

Aw = (A- X F(A(A - A I) + G) + A l)w . (27) 

The complex matrix eigenvalue problems above are solved using the LAPACK's 
ZGEEV routine Q. 
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— 38 947789 


— 66 057895 


— 60 054933 
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— 60 055435 


-73.670710 


-88.285123 


-88.299997 


-88.299997 




-119.43366 


-119.27480 


-119.27480 




— 157.89593 


— 157.38866 


— 157.38866 




-199.64318 


-198.23234 


-198.23234 




-226.99053 


-246.21576 


-246.21576 




-384.38914 


-296.92876 


-296.92874 




-409.06660 


-354.78191 


-354.78176 




-961.90740 


-415.36266 


-415.36420 




-961.99676 


-483.07721 


-483.08684 






-553.58796 


-553.53879 






-711.16464 


-711.45255 



Table 1: The eigenvalues found numerically by solving Eq. H21H (method I) with various 
number of collocation points N for a = 1 and Re = 0. The exact eigenvalues values are the 
roots of the characteristic equation resulting from analytical solution of Eq. {3} for Re = 0. 



5. Numerical results 

In order to validate the approach developed above we start with Re = 
for which Eq. Q can easily be solved analytically leading to the characteristic 
equation 

tanh(fc) =± ( k \ ±X / 2gl 
tantV^A) \Vk^~\J ' 1 ' 

which defines two branches of eigenvalues A for the even and odd modes corres- 
ponding to the plus and minus signs in the above expression. The eigenvalues 
resulting from Eq. (|2ip , which represents our first approach, are shown in Table 
[T]for various numbers of collocation points along with the exact solution defined 
by Eq. (|28p . As seen, this approach indeed produces a couple of huge spurious 
eigenvalues, which are due to the singularity of the transformation matrix E 
(|20p pointed out above. At the same time, the numerical solution accurately 
reproduces the leading eigenvalues of the exact solution. The accuracy, however, 
decreases down the spectrum so that only a half of the exact eigenvalues are 
reproduced by the numerical solution. The other half are numerically spurious 
eigenvalues which are due to the discretization of the problem Q • 

Our second approach defined by Eq. (f2~4"|) produces exactly the same eigen- 
values as the first one for the given N except for the two spurious eigenvalues 
which are now machine-size zeros rather than infinities. Using the modification 
of the second approach defined by Eq. ((27]) . these zero eigenvalues can be shifted 
to any prescribed value Ao without affecting other eigenvalues. Further we use 
Ao = 4(iV/4) 4 which shifts the two physically spurious eigenvalues to the region 
of numerically spurious eigenvalues located in the lower part of spectrum. The 



6 




Figure 1: Relative variation of leading eigenvalues with the number of collocation points N 
for a = 1, Re = and Re = 10 4 . 



N 


c 


16 


(0.23272286, 0.00922887) 


20 


(0.23814366, 0.00566303) 


24 


(0.23842504, 0.00282919) 


28 


(0.23735182, 0.00357013) 


32 


(0.23747200, 0.00372519) 


36 


(0.23752527, 0.00375400) 


40 


(0.23752494, 0.00373917) 


44 


(0.23752716, 0.00374012) 


48 


(0.23752633, 0.00373961) 


52 


(0.23752653, 0.00373969) 


56 


(0.23752648, 0.00373967) 


60 


(0.23752649, 0.00373967) 


64 


(0.23752649, 0.00373967) 



Table 2: Phase velocity c = — i\/(Rek) of the most unstable mode depending on the number 
of collocation points N for a = 1 and Re = 10 4 . 



variation of the five leading eigenvalues with the number of collocation points 
N plotted in Fig. [T]shows an exponential convergence rate characteristic for the 
spectral numerical methods [3j. 

Next, we consider the solution of Eq. (|27p for Re = 10 4 and a = 1, which is 
a standard test case for the linear stability analysis of plane Poiseuille flow. The 
leading eigenvalue for this case is shown in table [5] in terms of the commonly 
used phase velocity c = —iX/Rek. For N > 60 the solution is seen to converge 
to the reference value obtained in [l2| using a tau method with M > 30 even 
Chebyshev polynomials. As seen in Fig. [TJ however the convergence rate for 
Re = 10 4 is somewhat slower than for Re = 0, it is still exponential with the 
final accuracy comparable to the previous case. 

The number of collocation points can be reduced by a half by considering 
even and odd modes separately as done in [12j. In our case, this would require 
substitution of differentiation matrices (|29l30p for general functions with their 
half-size counterparts for even and odd functions, which, however, lies outside 
the scope of this paper. 
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6. Summary and conclusions 

We have developed a simple technique for avoiding physically spurious ei- 
genmodes in the solution of hydrodynamic stability problems by the Chebyshev 
collocation method, which was demonstrated on the Orr-Sommerfeld equation 
for plane Poiseuille flow. The method is based on the factorisation of the ori- 
ginal fourth order differential equation into two second-order equations using 
a vorticity-type auxiliary variable which has no explicit boundary conditions. 
The main element of the method is the elimination of the vorticity boundary 
values by using a capacitance matrix approach to obtain a standard matrix 
eigenvalue problem. Although the elimination is constrained by of the struc- 
ture of eigenvalue problem, it can be still done in two basically different ways. 
Both approaches result in couple of physically spurious eigenvalues, which are 
either huge or close to zero depending on the way the vorticity boundary val- 
ues are eliminated. We showed that these spurious eigenvalues are due to the 
double singularity of the transformation matrices which eliminate the vorticity 
boundary conditions by multiplying either the stiffness or mass matrices of the 
original generalised eigenvalue problem. By a slight modification of the second 
approach, the zero eigenvalues can be shifted to any prescribed value and thus 
moved to the region of numerically spurious eigenvalues at the end of spectrum. 

The main advantage of our method is not only its simplicity but also its 
applicability to more general stability problems with complex boundary condi- 
tions involving several variables. An example of such a problem is that of 3D 
linear stability of MHD duct flow using a non-standard vector stream function 
and vorticity formulation which results in the coupling of the stream function 
components through the boundary conditions [14]. In this case neither Galer- 
kin nor collocation method with the ad hoc approach of Huang and Sloan Q is 
applicable because no simple basis functions satisfying the boundary conditions 
can be constructed. 



Appendix 



The differentiation matrices for the first and second derivatives at Chebyshev- 
Gauss-Lobatto nodes Q are defined as follows [13], pp. 393-394]: 



and 



0$ = (A a J) a = < 



2N 2 + 1 
Qj (-1)' +J 

Ci (Zj-Zi) 



2(l-z 2 ) 
2N 2 + l 
6 



3 = 

< i = j < N 
j = N 



(_1)'+J z 2 +z i z j -2 
(AT 2 -l)(l-z 2 )+3 

m-z'y 

2 (-I) 3 (2jV a + l)(l- Zj )-6 

3 c 3 (1^7? 

2 (-iy + N (2N 2 + l)(l+z j )-6 

3 Cj (1+Zj) 2 
7V 4 -1 

15 



< % f j < N 

< i = j < N 
j^i = 
j^i = N 

1 = j = 0,N, 



(29) 



(30) 



where Ci = 1 for < i < N and ci — 2 for i = 0,N. 
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